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ABSTRACT 



We present a power spectral analysis of the Crab pulsar's timing noise, mainly 
using radio measurements from Jodrell Bank taken over the period 1982-1989, an 
CN \ interval bounded by sparse data sampling and a large glitch. The power spectral 

^ \ analysis is complicated by nonuniform data sampling and the presence of a steep 

^ \ red power spectrum that can distort power spectra measurement by causing 

00 ! severe power "leakage" . We develop a simple windowing method for computing 

^ \ red noise power spectra of uniformly sampled data sets and test it on Monte Carlo 

O ! generated sample realizations of red power-law noise. We generalize time-domain 

2 \ methods of generating power-law red noise with even integer spectral indices to 

"^H \ the case of noninteger spectral indices. The Jodrell Bank pulse phase residuals 

are dense and smooth enough that an interpolation onto a uniform time series 
O \ is possible. A windowed power spectrum is computed revealing a periodic or 

\ nearly periodic component with a period of 568 ± 10 days and a 1//^ power-law 

\ noise component in pulse phase with a noise strength S^f, = (1.24 ±0.067) x 10^^^ 

I cycles^/sec^ over the analysis frequency range / = 0.003 — 0.1 cycles/day. This 

^ ' result deviates from past analyses which characterized the pulse phase timing 

I residuals as either 1//^ power-law noise or a quasiperiodic process. The analysis 

was checked using the Deeter polynomial method of power spectrum estimation 
that was developed for the case of nonuniform sampling, but has lower spectral 
resolution. The timing noise is consistent with a torque noise spectrum rising 
with analysis frequency as / implying blue torque noise, a result not predicted 
by current models of pulsar timing noise. If the periodic or nearly periodic 
component is due to a binary companion, we find a mass function /(M) = 
(6.8 ± 2.4) X 10""^^ Mq and a companion mass, Mc > 3.2M0 assuming a Crab 
pulsar mass of 1.4M0. 
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1. Introduction 

Isolated rotation-powered pulsars generally exhibit a pulse frequency that can be mea- 
sured to high precision and that decreases, for the most part smoothly, on long time scales, 
indicating a steady loss of angular momentum from the rotating neutron star. The determin- 
istic spindown trend of the pulsars can be modeled by a low order polynomial in the pulse 
phase. Apart from well known "glitches" in rotation rate occasionally observed in a few 
pulsars, many pulsars show fluctuations in the pulse phase about the deterministic trend. 
These rotation fluctuations or "timing noise", often do not have the form of white noise 
but rather exhibit a "wandering" behavior with a roughly polynomial appearance. This is 
not a sign of unmodeled higher order polynomial terms in the deterministic spindown trend, 
however, since adding more polynomial terms does not improve one's ability to predict the 
pulse phase. The nature of these rotation fluctuations is of interest since they provide direct 
observational information on the forces acting upon the neutron star crust. The presence of 
timing noise also complicates the detection of perturbations to the pulse phase caused by or- 
biting neutron star companions and has the potential to be mistaken for such perturbations 
if not properly understood. 

The Crab pulsar spins at a rotation frequency of 30 Hz and is the most well observed 
of all isolated pulsars. The Crab pulsar was first discovered in 1968 and has been subjected 
to many long term observation programs at optical and radio wavelengths (Groth 1975a; 
Gullahorn et al. 1977; Lohsen 1981; Lyne, Pritchard & Smith 1988, 1993). There have been 
several prior analyses of the long term properties of the pulse phase of the Crab pulsar. They 
all agree that the intrinsic pulse phase can be divided into three components: (1) A long 
term spin down that can be characterized by a low order polynomial (2) several glitches and 
(3) timing noise. 

The nature of the Crab timing noise has been investigated a number of times. The 
timing noise is revealed as the dominant component in the timing residuals after removal 
of a spindown polynomial and glitches from the pulse phase. The timing residuals from 
optical observations were compared to noise models consisting of random walks in pulse 
phase, pulse frequency or pulse frequency derivative (Boynton et al. 1972; Groth 1975b, c; 
Cordes 1980) and found to be most consistent with a random walk in pulse frequency in all 
cases. A low resolution power spectrum was calculated for the pulse frequency derivative 
using a polynomial estimator method by Deeter (1981) (also displayed in Boynton 1981; 
Alpar, Nandkumar & Pines 1986), who found a white spectrum consistent with a random 
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walk in pulse frequency. For a review of this work see D'Alessandro (1997) and Deshpande et 
al. (1996). In contrast, using a five year span of radio observations from Jodrell Bank in the 
1980's, Lyne, Pritchard & Smith (1988) claimed that the timing residuals were consistent 
with a physically real quasiperiodic process with a period of ~ 20 months since the period 
did not scale with the length of the data being fit, as occurs with the apparent periods 
observed in the residuals of polynomial fits to red noise. The disagreement between these 
various analyses has motivated this study of the timing properties of the Crab pulsar and 
a search for analysis techniques to better characterize the statistical nature of the timing 
noise. 

We focus on power spectral analysis methods in this paper to characterize the Crab 
pulsar timing residuals. The power spectral analysis is not straightforward due to a steep 
rise in power at low analysis frequencies and due to the nonuniform sampling of the data sets. 
The Crab timing residuals most likely contain some type of red power-law noise process, of 
which the random walk is a well-known example. Random walks have power spectra that 
rise as 1/ P toward low analysis frequencies (/). We will refer to / as analysis frequency 
to prevent confusion with the pulse frequency v. In general, a steep rise in power at low 
analysis frequencies can cause severe power "leakage" into the low frequency side-lobes of the 
frequency response of power spectral estimators, such as the squared amplitudes of a Fourier 
transform, that can be comparable to or greater than the response at the nominal frequency 
of the estimator itself. The leaked power greatly distorts the computed power spectrum 
preventing detection of its true form unless special precautions are taken to control this 
leakage. 

The nonuniformity of sampling has been severe in the pulse timing data sets obtained 
from optical observation programs that were carried out during the interval from 1969 to 
1980. The optical observations are subject to seasonal and monthly interruptions due to 
conjunctions of the Crab pulsar with the sun and moon. The total length of a data set that 
can be analyzed using a single pulse phase ephemeris has been limited as well due to the 
occurence of occasional large glitches that hinder correct extrapolation of the pulse phase 
across the glitch. Long term radio monitoring of the Crab pulsar at Jodrell Bank from 1982 
to the present has produced a set of pulse arrival times with much smaller gaps than in the 
optical data and generally has a higher sampling rate. In this paper we concentrate on a 
subset of the Jodrell Bank data set over the 7.5 year interval from February 1982 to August 
1989 that is bounded by sparse data samphng prior to 1982 and a large ghtch that occurred 
in 1989. 

Mathematical characterization of the timing noise in the Crab and other pulsars has 
frequently been in terms of red power-law noise processes that consist of random walks 
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in the pulse phase, pulse frequency or pulse frequency derivative or mixtures thereof (e.g 
Cordes and Downs 1985). One of the major attractions of these models is that discrete 
random walk realizations (and summations thereof) can be easily generated numerically 
using time domain methods. Random walks and their repeated summations (or integrals in 
the case of continuous sampling) all have ensemble averaged power spectra of the form 1//™ 
where the power-law spectral index m is an even integer. These noise models have been 
important for checking and calibrating various red noise analysis methods. We show that 
these time domain methods can be generalized to produce red-power law noise processes 
with arbitrary spectral indices (see appendix A). In addition, red (and blue) power law noise 
realizations with arbitrary spectral indices can also be generated using frequency domain 
methods. We use both time domain and frequency domain simulated noise realizations to 
check the accuracy, efficacy and calibration of the analysis procedures used in this study. 

Analysis of pulse arrival times necessarily requires the fitting of a low order polynomial 
to obtain timing residuals. When the timing noise has a steep red power spectrum, the 
polynomial fitting removes noise power at the lowest frequencies, creating a low frequency 
turnover or peak in the power spectrum. In the time domain, this tends to produce residuals 
with a quasiperiodic appearance. The appearance of a low frequency quasiperiodic looking 
residue can easily be mistaken for a real quasiperiodicity or even a real periodicity by the 
unwary investigator when in fact the quasiperiodicity is purely an artifact of the polynomial 
fitting and the quasiperiod will scale with the length of the data fit. To detect a true 
low-frequency quasiperiodicity in the presence of detrended red noise requires more careful 
analysis than simple visual inspection of the timing residuals. We show that in the Crab 
pulsar, a real periodicity or quasiperiodicity is present in the power spectrum of the Crab 
timing residuals produced from the Jodrell Bank observations in addition to a red noise 
component. 

The analysis of Deeter (1981) involved developing and applying a method of computing 
low frequency resolution power spectra that deals with the problems of nonuniform sampling 
and red noise leakage. The method uses polynomial fits on a hierarchy of time scales to 
generate a power spectrum at octave frequency resolution. The low frequency resolution 
of this method means that any significant narrow features in the power spectrum will be 
unresolved and allows significant room for misinterpretation of the power spectral shape. 
We show that high resolution power spectra can be computed for the case of uniform data 
sampling by use of a simple windowing procedure in conjunction with polynomial detrending 
to alleviate the leakage problem. The smoothness and dense samphng of the Crab timing 
residuals produced from the Jodrell Bank observations makes possible an interpolation onto 
a uniform grid and computation of a high resolution power spectrum. This new analysis 
shows that the pulse phase timing noise is composed of two components: a red noise process 
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with an approximately 1/P power density spectrum and a long term periodic or nearly 
periodic oscillation with a period of 568 ± 10 days. We confirm this result, at lower spectral 
resolution, using the Deeter polynomial method. 

2. The major components of the pulse phase of the Crab pulsar 

The pulse phase of the Crab pulsar is obtained via the application of a pulse phase 
ephemeris to a set of measured pulse arrival times that have been corrected to the solar 

system baryccntcr and corrected for other timing delays such as those caused by interstellar 
dispersion of radio waves (see Lyne, Pritchard & Smith, 1988). For the Crab pulsar, the 
intrinsic pulse phase can be decomposed conceptually into three components: 1) the pulsar 
spin frequency and its deterministic decrease, 2) glitches and 3) rotation fluctuations about 
the deterministic trend. These appear to be essentially independent components although 
the presence of glitches may be effecting the long term slow-down process (Lyne, Pritchard & 
Smith, 1993). In addition, we will consider the possible presence of a longterm quasiperiodic 
or periodic variation in the pulse phase. The pulse phase can be represented by a model: 

(j){t) = (l)s{t) + (j)G{t) + (j)TN{t) + (j)Qp{t) + (j)M{t) (1) 

representing respectively the contributions of the rotation frequency and its deterministic 
spindown (^5), glitches (</>«)) intrinsic timing noise (^tat), a quasiperiodic process (0qp) 
and measurement errors (0m)- In order to study the properties or detect the presence of 
any random or longterm periodic processes in the pulse phase we have to investigate the 
properties of the pulse phase residuals 4>res{t) = 4>{t) — 4>eph{t) on a time span T rather than 
(j){t) directly. The pulse phase ephemeris usually consists only of a low order polynomial to 
model 05, although "glitch functions" are sometimes added to fit the ghtches. For the data 
sets analyzed in this paper, glitches will either be absent or only a minor perturbation on 
the pulse phase residuals that remains after subtraction of the pulse phase ephemeris. 

The deterministic spindown of the Crab pulsar is modeled using a simple spindown 
model, z> = —Ki/^, where u is the rotation frequency and n is the braking index. The value 
of the braking index depends on the mechanism slowing down the pulsar. For braking purely 
by the emission of magnetic dipole radiation n — 3 and for braking only by emission of a 
wind n = 1. The braking index can be determined using n — if i>o can be measured. 
The analytic solution to the spindown model is given by: 

(2-n) 

Mt - to) = i^o'^oC^ - n)-' ([1 + ^>oi^o"'(l -n){t- to)] ^ - 1) + 0o. (2) 
The term uoi^q^ is quite small for all known pulsars so the analytic solution can be expanded 
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using the binomial theorem into a "spindown polynomial" : 

(Ps{t - to) = 00 + Mt - to) + - to)^ + - to)^ + ^ ^0 - to)^ .... (3) 

The Crab pulsar has a measured pulse frequency derivative of z>o = —3.862 x 10~^° Hz/s and 
a braking index of 2.5 (Lyne, Pritchard & Smith 1993) thus successive coefficients in the Crab 
spindown polynomial will decrease in magnitude by a factor of 10~^°. This implies that, if 
the rotation phase consisted only of this spindown polynomial, a polynomial fit of order I will 
always produce polynomial residuals of order / + 1, unless masked by measurement errors. 
The timing noise in the Crab pulsar was initially revealed by the presence of polynomial-like 
residuals that were much larger those expected from a pure spindown polynomial (Boynton 
et al. 1972). 



3. Observations and computation of pulse timing residuals 

A set of radio pulse arrival times for the Crab pulsar corrected to the solar-system 
barycenter were obtained from Jodrell Bank (Lyne, private communication). The data con- 
sist of 7177 arrival times collected over a nine year span from January 27, 1982 to January 
23, 1991 (Truncated Juhan Days 4996 to 8279; TJD = JD - 2440000.5) at five radio frequen- 
cies (408, 610, 930, 1420 and 1667 Mhz) with approximately 85% of the arrival times at 610 
Mhz. The typical arrival time error was 9 microseconds. These data or portions thereof have 
been previously presented in Lyne and Pritchard (1987), Lyne, Pritchard & Smith (1988) 
and Lyne, Smith & Pritchard (1992) and Lyne, Pritchard & Smith (1993). Two glitches are 
present in the data occuring on TJD 6664 with a Au/i/ ~ 10~^° (Lyne and Pritchard 1987) 
and TJD 7767 (Lyne, Smith and Pritchard 1992). The second glitch is one of the largest 
reported for the Crab pulsar {Au/u ~ 10~^) and the post glitch arrival times are sufficiently 
distorted from their pre-glitch behavior that only data prior to this glitch were included in 
the subsequent analysis (4530 arrival times). In contrast, the glitch occuring on TJD 6664 
was a relatively small event and is discussed in section 8. 

A quartic polynomial phase model of the form given by equation 3 was removed from 
the arrival times using values given in tabic 1 and the method described on page 104 of 
Manchester and Taylor (1977). The errors in the parameters were determined by the method 
of fitting a quartic polynomial -|- sinusoid to the phase residuals determined from an initial 
phase ephemeris. See section 8 for a description of this process. 

The power-law spindown model given by equation 2 predicts a value for z/q of 

i>o = n(2n - l)z>o/i/o 
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= -6.2 X 10"^^ Hz s-^ (4) 

which is approximately half the value in table 1. The phase contribution of the measured uq 
term is 129 cycles over a period of 2500 days, whereas the value predicted by the power-law 
spindown model would contribute 56 cycles. Hence z/q should be easily measureable without 
much distortion caused by red noise interference. The secular spin-down thus appears to be 
somewhat more complex than that predicted by a simple power-law spin-down model. This 
may be an indication of the longterm effects on the spin down process caused by glitches 
(Lyne, Pritchard & Smith (1993)). 

The arrival time data recieved from Jodrell Bank had been corrected only for a single 
value of the dispersion measure (DM) over the entire data span. Therefore, prior to removal 
of the quartic polynomial phase model, the arrival times were corrected using the differ- 
ence between monthly DM measurements taken from the Crab pulsar monthly ephemeris 
available at the world wide website: http:/ /www .jb.man.ac.uk/ pulsar/crab. html and a sin- 
gle DM value of 56.743. The changes in the DM value from month-to- month caused small 
discontinuities in the arrival times at the bin boundaries. This problem was dealt with by 
making a cubic sphne interpolation to the set of monthly DM values at a higher time reso- 
lution and then computing DM values with a linear interpolation at the actual data times 
using the nearest pair of DM values from the cubic spline interpolation. The resulting timing 
residuals are quite smooth and very similar to that displayed in figure 6 of Lyne, Pritchard 
& Smith (1993). A fair number of points (130) deviated significantly from the local trend 
in the timing residuals. These were judged to be outliers and removed. We show the timing 
residuals in the top panel of figure 1. The timing residuals arc relatively smooth and dis- 
play an apparent quasipcriodic process. The number of zero-crossings is ten, five more than 
expected if a pure spindown polynomial were present. 



4. Red Power-law noise processes 

The timing noise in the Crab pulsar has been attributed to a random walk in pulse 
frequency by several previous analyses (Boynton ct al. 1972; Croth 1975b, c; Cordes 1980). 
A random walk is a member of the class of power-law noise processes and has an 1//^ 
power-law power spectrum. In pulse phase, that is the integrated pulse frequency, the power 
spectrum steepens to a 1//"^ power-law. We will refer to these type of noise models as 
red power- law noise (PLN). The a priori knowledge that red PLN is probably present is 
important since analysis methods and interpretations that usually assume only white noise 
is present would produce erroneous results. 
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Random walks in pulse phase, pulse frequency and pulse frequency derivative^ have been 
considerably investigated with respect to pulsar timing noise since these are simple models 
with which noise realizations can easily be generated. They are members of the class of red 
PLN processes with 1//™ spectra where m = 2,4, 6, . . . that can be generated by repeated 
integration of white noise. Integration of white noise processes will produce random walks 
which have an ensemble averaged power density spectrum with a power-law spectral index 
of m = 2. Integration of the random walks will produce noise realizations with an ensembled 
averaged power-law spectral index of m = 4. The time series becomes smoother with each 
integration and hence "redder" in appearance as the low frequency components of the time 
series become more dominant. For example, a white noise process in the pulse frequency 
derivative will produce a random walk in pulse frequency and an integrated random walk 
in pulse phase. In general, each integration will increase the power-law index of the power 
density spectrum of the noise process by 2. To see this, consider a zero mean continuous^ 
white noise process, e(t), with an autocovariancc function < e(t)e(t + r) >= a^5(r), where 
6{t) is the Dirac delta- function. The quantity S = a"^ is called the "noise strength" and is 
a fixed parameter that characterizes the white noise. The ensemble averaged power density 
per unit analysis frequency of e(i) will be a constant, < \F^{f)\'^ >— S where / is analysis 
frequency and F^{f) is the Fourier transform of e{t). The derivative theorem for Fourier 
transforms states that: if f{t) has a Fourier transform F(f) then f'{t) has a Fourier transform 
i27r/F(/). Let us consider a random walk defined by: 

r,{t) = r e{t)dt (5) 

where the subscript on r{t) will refer to the power-law spectral index of the process which 
has Si. 1/ p ensemble- averaged spectrum and the Fourier transform of r2(t) is Fr^if)- Since 
€{t) is a continuous random process existing at all t on — oo < t < +oo, and F^(f) exists 
in the limit of the Fourier Transform of band-hmited white noise, the fact that r'2{t) — e(t) 
imphes that the Fourier transform of r'2 is i2T:fFr^{f) and is equal to F^{f). Therefore, 
Fr^if) = Fe(/)/i27r/. The ensemble averaged power density spectrum of r2{t) is defined 
by: Pr-zif) = -Pe(/)/47r^/^ = S/An'^p. This argument can be repeated for each successive 
integration. Thus each integration will steepen the power density spectrum by a factor 
l/(47rV'). 



^Groth (1975b) refers to these respectively as "phase noise" (PN), "frequency noise" (FN) and "slowing- 
down noise" (SN), whereas Deeter (1984) refers to 1//™ power-law noise with m = 2,4, 6 as 1st, 2nd and 
3rd order red noise respectively. 

^Continuous in the random variable sense, i.e. defined for all times t 
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Previous timing analysis of the Crab pulsar (Boynton et al. 1972; Groth 1975b, c; Cordes 
1980 and Deeter 1981) focused on power-law noise processess that could be generated by 
the integration or summation of white noise realizations and so was biased toward power- 
law noise processes with even integer spectral indices. There is no a priori reason that 
the timing noise should be restricted to an even integer red power-law noise process and 
a truly general analysis method should be able to tell the spectral index of any power-law 
noise directly rather than making a best fit from among a set of expected indices as has 
been done in several of these past analyses. The Deeter polynomial method (Deeter & 
Boynton 1982, Deeter 1984) can in principle determine the power-law index directly but 
has such low spectral resolution that the presence of any discrete features such as those 
due to a quasiperiodic process will produce an erroneous interpretation. A need exists for 
higher spectral resolution methods that can directly determine the power spectral index and 
distinguish the presence of narrow features, which we explore below. 

In order to test possible analysis methods, we develop a time domain method to generate 
red power-law noise realizations with an arbitrary spectral index in Appendix A, which is 
a generalization of the discussion given by Groth (1975b). A useful, but unnormalized, 
frequency domain method which weights Fourier amplitudes by a power-law and then inverse 
Fourier transforms the result into the time domain has been presented by Timmer & Konig 
(1995) which also generates PLN realizations with an arbitrary spectral index and is very 
straightforward to implement (see Appendix B for a normalized method). 



5. Power Spectral Density Estimation of Red Noise Processes 

The estimation of the power density spectrum (PDS) of a red PLN realization presents 
special problems due to the steep power rise in the spectrum at low frequencies. The problem 
arises due to leakage of power in the low-frequency sidelobe of the frequency response of a 
spectral density estimator which can greatly exceed the response at the central or nominal 
frequency. Consider a noise realization r^(t) over the time interval ^ to The estimated 
PDS will be the squared modulus of the convolution the Fourier transform of the "window 
function", W{f) with the intrinsic Fourier transform of the noise process F{f): 

/+00 
F{f)W{f - f)df (6) 
-oo 

If the time interval extended over infinity then the Fourier transform of the window function 
would be a 5-function and the intrinsic power spectrum would be recovered exactly. The 
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ordinary rectangular window function: 
has a Fourier transform given by: 



1 for \t\ < T/2 
otherwise 



When F{f) has a power- law form: (27^y^, the windowed Fourier transform is given by: 

C sin Wif - f')T] , 

The power leakage problem occurs due to the singularity in the power-law noise spectrum 
at zero analysis frequency. 

In practice, the rise in power towards zero frequency will exhibit a cut-off due to the use 
of a finite, rather than infinite, stretch of data. The main leakage problem is caused by the 
large, but finite power, at the lowest frequencies actually present. When m > 2, i.e. for red 
PLN as red or more red than a random walk, each frequency estimate will be dominated by 
the response to the power at the lowest frequencies. The estimated spectrum will then have 
a power-law spectral index of ~ 2 regardless of the value of m. Detrending the data prior to 
power spectral estimation, by a polynomial fit for example, will decrease the power at the 
lowest frequencies and will help to alleviate the leakage problem irregardless of the type of 
power spectral estimator applied to the data with the tradeoff that the nondetrendcd form 
of the power spectrum cannot be measured directly. These issues will be further explored in 
section 7. 



6. Computation of a Windowed PDS in the Presence of Red Power-Law Noise 

A PDS of a uniformly sampled red noise time- series can be computed with fast Fourier 
Transforms by using an appropriate window function to control the power leakage. The 
window function must be chosen to sufficiently reduce the side-lobe response of the power 
density estimator (at the expense of spectral resolution) so that leakage is not a major 
problem. For a given time-series, {xk, tfc} A; = 0, 1, ... — 1, the Hann window of the form 
Wk — sin°'{7Tk/N) with a — 2 has a frequency response with sidelobes whose envelope falls- 
off with analysis frequency as 1//^ (see Harris 1985). This is sufficient for computation 
of power density spectra for red noise with power-law indices up to m ~ 5 according to the 
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moment condition criteria (Deeter & Boynton 1984). Presumeably a window function could 
be designed to optimally suppress the side-lobe response while preserving maximum spectral 
resolution but this has not been pursued in this paper, since sine windows have proven to 
be simple and adequate in practice. 

The actual method of PDS computation and normalization used here involves several 
steps. After a window, Wk, is chosen the window is normalized to a root-mean square of 1: 



Wk 



The PDS of WkXk is computed as: 



N-\ 

WkXkexp{+2nijk/N) 

k=0 



. _ N N 



By taking advantage of the fact that for real data \ aj\ — \ a-j \ , the two-sided PDS is converted 
to a one-sided PDS, ja^p, by doubhng all powers except |aoP and the Nyquist term (|a_ivp) 
when N is even and eliminating all redundant negative frequency terms. 

A normalization factor of is applied to the \aj\^. This is done so that if Xk is a white 
noise process then the mean power level will equal the white noise variance. For example, if 
Xk is a white noise process with < x\ >— cr^ and < Xk >— then: 

AT-l 




k=0 
N-1 

W 



k 
k=0 
2 



Nat (8) 

(9) 



Using Parseval's relation with the constant term (j = 0) deleted we have: 

1 



N 



/ (10) 

(11) 
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where \aj\'^ are the one-sided powers defined above and P is the mean power level. We 
therefore define a normalized PDS as: 

P = —laA^ j = l — (12) 

so the mean power level will approximately equal the white noise variance cr^ if Xk is a 
zero-mean white noise variable. 

The power-law index of the red noise noise process can be determined by performing a 
linear fit to log Pj vs. log /_,■ where the fj = are the analysis frequencies. An appropriate 
frequency range must be chosen for the fit to avoid frequency regions where the power 
spectrum contains substantial contributions from any other processes present in addition to 
a pure power-law noise process. In this paper we used linear least-squares fitting to estimate 
the slope. The data points log Pj do not have a gaussian distribution about a mean value 
< log Pj > as assumed in ordinary least-squares fitting, so some distortion in the slope and 
slope error estimates will be introduced. The actual distribution of the data points about a 
mean value is similar to gaussian but possesses a low value tail that results from taking the 
log of a chi-squared distribution with two degrees of freedom. The actual slope estimation 
procedure used here involved an initial fit with assumed unit errors on log P, a calculation 
of the deviation of the fit residuals and then a new fit with uniform errors equal to the 
measured deviation. Monte Carlo simulations of this procedure on simulated data sets of 
10000 sample realizations each for a variety of input power-law indices showed that no bias 
existed in the estimated mean slope value but the slope errors had been overestimated by a 
factor of 1.7. Hence, the slope errors quoted in this paper have been divided by 1.7. 

The power density spectrum can be converted to a flattened form to allow clearer study 
of features unrelated to the power-law form. We chose a normalization that allows the 
noise strength to be estimated directly from the mean power-level over a suitable analysis 
frequency range. Given a power-law index m the power spectrum is multiplied by {2'n'fj)"^. 
The analysis frequency is converted to physical units as = (^Ai^ where At = — 
The flattened PDS is then normalized so the mean power level approximately equals the 
noise strength in physical units by multiplying by a factor {-^)"^~^- This is a combination 
of a factor of (;^)"* to flatten the spectrum, a factor of (;^)~^ to correct the power via the 
similarity theorem for the rescaling of the frequency and a factor of (^) to make the mean 
power level equal the noise strength in the presence of pure power-law noise. The extra 
factor of (■^) arises from the deflnition of the noise strength S (see section 7 and Appendix 
A). The total set of operations on the one-sided PDS can be combined as: 



j, flattened 



1 \ / 1 \"'"', .„ . _ iV 
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Note that the power-law spectral index m is a free parameter supplied by the user. The 
units of the noise strength are given by the units of Uj and At. For example, a random walk 
in pulse phase (measured in pulse cycles with time measured in seconds) will have m = 2 
and hence a noise strength S = {P flattened) with the physical units of cycles^ /s, whereas a 
random walk in pulse frequency (m=4) will have the physical units of cycles^ /s^. 

7. Power Spectrum Estimation of Simulated Red Noise Realizations with 

Uniform Sampling 

We explored the practical estimation of power spectra of red PLN processes using sim- 
ulated PLN time-series. This problem involved determining methods to generate PLN pro- 
cesses and then seeing if the predicted spectra and noise strength could be correctly recov- 
ered. Neither the generation of red PLN processes with non-even integer spectral indices 
nor high resolution estimation of their power spectra has been previously well developed, 
so sorting out the presence of possible systematic effects and attributing them to either the 
noise generation method or the spectral estimation method required extensive testing. 

We implemented the time domain method of generating discrete red PLN realizations 
described in Appendix A. With this method a band limited white noise realization e{t) 
is used as a basis upon which a power-law weight function acts to generate the red noise 
realization. The white noise is a zero mean process with a gaussian distribution of impulse 
amplitudes generated at times tk with nonuniform time steps (At — tk+i — tk), where At 
is drawn from the exponential waiting time distribution of a Poisson process. The chosen 
step rate Rg and variance of the step amplitude distribution cr^ determine the noise strength 
5* = RgC^- The power-law weight function then acts upon the discrete white noise realization 
to produce a set of accumulated steps that can be sampled at an arbitrary set of times ti that 
are independent of the times of the white noise impulses. For the case of a random walk, 
one will have a set of rectangular steps of varying amplitude and duration, whereas for 1//^ 
noise, the individual steps will have a form oc ^2. The appearance of the simulated noise 
realization for the same noise strength and power-law index varies appreciably depending on 
the sampling and step rates. For example, sampling at a rate greater than the white noise 
step rate Rg allows individual steps to be resolved. In figure 2 we portray a set of discrete 
red power law noise samples with resolved steps. 

In figure 3 we show a set of red noise realizations generated from the same underlying 
white noise realization. As the power-spectrum becomes "redder" the time series takes on 
an increasingly smooth appearance that arises from the increasing dominance of the lowest 
frequency components. The right-hand set of panels shows the same noise realizations after 
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cubic polynomial detrending. A quasiperiodicity is apparent in the residuals in some cases 
that is a result of the polynomial detrending. Polynomial fits on shorter or longer timescales 
produce similar looking residuals and hence a "pseudoperiodicity" that scales with the length 
of the red noise time-series. 

In the top panel of figure 4 we show an example of a 1//^ noise reahzation that has been 
detrended by a cubic polynomial fit. All the time-series discussed in this section are uniformly 
sampled at a rate approximately twice the step rate and the results discussed apply strictly 
only for uniform sampling. In the middle panel we show the computed power spectrum 
(rough histogram) computed using a Hann window with a = 2. Lying above this is a 
smoother histogram that shows the average power spectrum of 1024 similar noise realizations 
(shifted in power by 3 decades for display purposes). The simulations showed that the input 
power spectral slope and noise strength could be correctly recovered. The spectrum shows a 
clear turnover at low frequencies caused by the removal of low frequency power by the cubic 
polynomial detrending and this results in residuals that appear quasiperiodic. The average 
analysis frequency of the peak is purely related to the length of the data span fit. In the 
bottom panel we show the same power spectra after having been flattened as described in 
section 6 assuming a power-law index m = 3. We now see that the effect of the polynomial 
detrending is to decrease low frequency power in the flattened spectrum relative to higher 
frequencies. A high frequency upturn in the average spectrum is also present that is the 
result of aliased power from above the Nyquist frequency. 

Examination of indvidual power spectra showed that cubic detrended realizations with 
a low number of zero-crossings (4-5) often displayed a slight excess of low frequency power 
with respect to the higher frequencies whereas those with higher numbers of zero-crossings 
tended to show a power deficit. This is an easily understood byproduct of the goodness 
of the cubic polynomial fit to the noise realization (i.e. better fits cause larger numbers of 
zero-crossings in the residuals). A second related effect is to produce a negative correlation 
of the mean squared value of the residuals with the number of zero-crossings, confirmed 
by simulations and is a general result for polynomial detrending of red noise. The degree 
of the polynomial fitting exhibited a weak positive correlation with the frequency of the 
power peak in the unflattcned spectrum. This is expected since quartic and higher order 
polynomial detrending simply removes more low frequency power and hence shifts the power 
peak frequency to higher frequencies and results in a greater deficit of low frequency power 
in the flattened power spectrum. 

Several effects of windowing and polynomial detrending are illustrated in figure 5 which 
displays averages of 1024 flattened power spectra for PLN realizations of unit strength with 
added white noise and a high frequency sinusoid. The power spectra were made using a Hann 
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window with a — 2. First, the polynomial detrended power spectra show clearly the decrease 
in low frequency power relative to the power spectra produced from undetrended noise 
realizations. The relative power drop increases as the power-law index increases consistent 
with the fact that redder realizations are smoother and thus better fit by polynomials. 
Secondly, the undetrended power spectra show an excess of power at low frequencies which 
increases with the power-law index. This is probably a power leakage effect caused by leakage 
into the highest sidelobes of the power estimator which will be largest at low frequencies 
where most power is concentrated. Thirdly, another leakage effect is also readily apparent 
in figure 5 illustrated by the spectra of the undetrended realizations: namely that for power 
spectral estimation at higher frequencies power leaks into the entire low frequency side of 
the estimator and hence tends to get worse as the estimator's nominal frequency increases. 
This effect eventually causes the spectral slope to change from the predicted power-law to a 
much lower value 2) when the increasing steepness of the power density spectrum causes 
the leaked power to exceed the rate of sidclobe fall-off of the power spectral estimator which 
for the Hann window with [a = 2) is ~ 1//^. The polynomial detrending, by decreasing the 
lowest frequency power, which is the largest source of power leakage, decreases substantially 
both the local sidelobe leakage and the integrated leakage on the low frequency side of the 
estimator. 

The high frequency area of figure 5 illustrates the dominance of the added white noise 
component at high frequencies as the red noise power level drops below that of the mean 
white noise power level. The "crossover" frequency is that power level where the mean red 
noise power level equals that of the white noise and occurs in this case near the frequency 
/ = 0.2. The high frequency sinusoid added to the noise realizations clearly shows up as 
the narrow peak near / = 0.1. In general, the correct sinusoid frequency could be recovered 
from the flattened spectra with no more bias than in ordinary power spectral estimation. 

We also tried other windowing functions such as the Hann window with a — 1. This 
window has better frequency resolution but worse leakage protection than for a — 2. Major 
leakage effects began to be noticeable for power-law indices for m > 3.5 for power spectra of 
nondetrended noise realizations and for m > 4.4 for detrended realizations. This is expected 
since the sidelobes fall-off as ~ 1//"^- Other window functions such as the Hamming, Kaiser- 
Bessel or Bohman were examined as well but rejected due to their poor leakage protection. 
The Blackman-Harris window offered leakage protection comparable or better than the cosine 
windows in some frequency regimes but had noticeably worse frequency resolution due to 
the large width of the main lobe of the frequency response and so was rejected as well. 

A second method of generating red PLN has been proposed by Timmer & Konig (1995) 
using a method that creates PLN realizations by generating the real and imaginary values of a 
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set of Fourier amplitudes and then weighting them by a power-law. The amplitudes are then 
Fourier transformed back into the time-domain. The use of a second method of generating 
PLN was useful in helping to differentiate effects caused by the spectral estimation technique 
versus those caused by the noise generation method. Direct comparison between the two 
noise generation methods was at first complicated by the lack of a proper normalization 
for the Timmer & Konig generated realizations. We determined a normalization for this 
method to generate noise with a predetermined noise strength consistent with that used for 
the time-domain method. See Appendix B for a more complete description. 

The Timmcr-Konig method generates a timc-scrics from a discrete set of Fourier fre- 
quencies and is therefore missing power from frequencies not within the discrete set. Because 
of this, there is no power at frequencies higher than the Nyquist frequency, or lower than 
the lowest frequency 1/T, where T is the total data timespan, as there would be in a more 
reahstic segment of PLN. If the time-series is resampled at a higher rate, then a small ahas- 
ing effect can be observed in the power spectrum at the highest frequencies. Likewise, if the 
time-series is resampled at a lower rate then a deficit of power is apparent at the highest 
frequencies. If one uses the full time-series, then the endpoints and the derivatives of the 
lowest order frequency terms nearly match when the time series is cyclically repeated, which 
greatly reduces the power leakage problem. Thus we discarded the second half of the time- 
series to better represent true power-law noise. Without discarding the second half, flattened 
power-spectrum with the correct power and slope can be produced using only a rectangular 
window for power-law indices much greater than m = 2, e.g. upto at least m — 5.5! This 
is just a result of artiflcially constructing red noise only from harmonics that are exactly 
periodic on the analysis interval which causes all leaked power to fall exactly on the zero's 
of the sine function frequency response. 

Average power spectra similar to that of figure 5 were calculated for noise realizations 
generated by the Timmer-Konig method. These power spectra, similarly to those in figure 
5, showed an excess in low frequency power for undetrended realizations and a subsequent 
drop in the power as a result of dctrending. Major leakage effects set in at similar power-law 
indices. The only real differences were the lack of aliased high frequency power and different 
high frequency distortion effects for power-law indices m > 5.5 where the Hann window with 
a — 2 ceases to provide adequate leakage protection. 
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8. Analysis of the 1982-1989 Crab radio timing residuals 

8.1. Windowed Power Spectrum 

We calculated a windowed power spectrum of the Jodrell Bank Crab timing residuals. 
The dense sampling and relative smoothness of the timing residual suggested that interpo- 
lation onto a uniform grid was possible without significant distortions. The timing residuals 

were first averaged on a one-day time scale, mainly to reduce the number of points in the 
densely sampled region after the 1986 glitch (see Lyne, Pritchard & Smith 1993 and figure 
6, middle panel). The data were then interpolated onto a uniform grid using a cubic spline 
interpolation with At = 1.346 days. Two samples comparing the uniformly gridded phase 
residuals with the unprocessed phase residuals are shown in figure 7. Using the method 
described in section 6 a power spectrum was calculated using a Hann window with a — 1. 
The uniformly binned time-series, windowed power spectrum and flattened power spectrum 
are displayed in figure 6. The power spectrum of the timing noise shows a power-law with a 
spectral index of m = 3.09 ±0.05. The flattened power spectrum shows a clear low frequency 
peak at a frequency / ~ 0.0018 cycles/day. In addition, the flattened power spectrum ex- 
hibits a decrease in low frequency power due to the removal of the polyonomial trend and an 
upturn at high frequencies due to white noise in the pulse phase. A power spectrum was also 
calculated using a Hann window with a = 2 (see bottom panel of figure 1) which provides 
better protection against red noise leakage but worse frequency resolution. Essentially the 
same power spectrum was produced but with worse resolution of the low frequency peak. 

The power spectrum of the Crab timing residuals is consistent with the presence of a 
power-law noise component with a spectral index m ~ 3. The noise strength is given by 
the mean white noise level in the fiattened spectrum. The averaged power in the fiattened 
spectrum over the frequency range 0.003 — 0.1 cycles/day is (9.26±0.5) x 10~^ cycles^/day^. 
This imphes a noise strength S — (1.24 ± 0.07) x 10~^^ cycles^/sec^. The timing noise is 
consistent with a torque noise spectrum rising with analysis frequency as ~ /, implying 
blue torque noise, rather than white torque noise as found by Deeter (1981), Cordes (1980) 
and Groth (1975c). For the numerical 1st derivative of the uniformly binned time-series 
we repeated the power spectral analysis described above and found a spectral index m = 
1.12 ± 0.05 and a mean flattened noise power (assuming m = 1) of (8.96 ± 0.5) x 10~^ 
Hz^. Likewise, the numerical 2nd derivative had a spectral index m — —0.86 ± 0.05 with a 
mean fiattened noise power (assuming m — —1) of (8.77 ± 0.5) x 10~^ (Hz/s)^day^. Both 
these results are consistent with that for the pulse phase. A large peak in the power showed 
up clearly at the same analysis frequency in the fiattened spectra for all three cases. The 
observed power spectrum of the pulse phase residuals is consistent with the presence of a 
power-law noise process, but we point out that this does not mean a power-law noise process 
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such as those described in appendix A is the correct form for the noise since there may be 
other noise processes that could produce similar power spectra. 

The cubic spline interpolation onto a uniform grid could introduce some distortions in 
the computed power-spectrum, in spite of the relative smoothness of the Crab timing residu- 
als. We checked the correctness of the analysis procedure described above using Monte Carlo 
simulations. A noise sample was generated by using the times of the Crab timing residuals 
after outliers had been removed. Using these sample times we generated a realization of 
1/ f"' power-law noise using the procedure described in Appendix A with a noise strength 
5* = 1.0 X 10^^ cyclcs^/day^ and a step rate Rg = 0.7434 day^-*^. To this we added a white 
noise realization with a variance equal to the variance of the Crab timing measurement noise 
and a cosinusoid with parameters given in table 2. The realization was then fit by a quar- 
tic polynomial, averaged on a one day timescale and interpolated onto a uniformly gridded 
time-series using a cubic spline following the same procedure used in analysis of the actual 
Crab timing data. Power spectra of the resulting time-series were computed using Hann 
windows with a — 2. 

We generated a sample run of 512 simulated noise realizations, computed power spec- 
tra, and measured the power spectrum slope and noise power level in the flattened power 
spectrum over the frequency range 0.003 to 0.1 cycles/day. This frequency range avoided 
the low frequency power contribution of the cosinusoid and the high frequency white noise 
contribution. We found a mean PDS slope of —3.03 ± 0.15 and a mean power level of 
(0.97 ± 0.1) X 10~^ cycles^/day^, consistent with the input values. The peak frequency of 
the input cosinusoid was correctly recovered. A second similar run of 512 sample noise real- 
izations without the low frequency cosinusoid was created and analyzed. A mean PDS slope 
of —3.03 ± 0.14 and a mean power level of (0.98 ± 0.1) x 10~^ cycles^/day^ was found. Any 
distortions due to aliasing would show up in the high frequency portion of the spectrum 
where the effects of white measurement noise tend to dominate and hence this is not a useful 
portion of the spectrum in any case. We therefore conclude that no significant distortions 
to the power spectra were introduced by the data analysis procedure described above. 

8.2. Quasiperiodicity 

The large low frequency peak in the Crab power spectrum appears to be a real feature 
and not simply a result of polynomial detrending of red noise. Detrended red noise should 
show a significant decrease in the low frequency power in the flattened spectrum relative to 
higher frequencies rather than an increase, as figure 4 demonstrates. The quartic polynomial 
removed from the Crab timing measurements will cause an even greater drop in low frequency 
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power in the flattened spectrum relative to that shown in figure 4 where cubic detrending was 
used. To investigate the significance of the power spectral peak we assumed the "white noise" 
in the flattened power spectrum followed a chi-square distribution. To check this assumption 
we made a histogram of the flattened power over the frequency range / = 0.003 — 0.1. The 
powers did indeed follow an exponential distribution. The peak at / = 0.0018 cycles/day 
has a probability of chance occurence of less than 1 x 10~^. 

To determine whether the power spectrum peak was due to a periodic or quasiperiodic 
process we performed a sinusoidal flt to the Crab phase residuals in the low frequency region 
around the power spectrum peak. We represented the Crab phase residuals in the time 
domain with a model given by 



•^model 



-^jp—^ + Acos{27rfAt) + Bsm{27rfAt) (14) 



where / is the sinusoid frequency and At — t — tg where to ~ TJD 6390 is an epoch. This 
model was fit to the data by doing a modified-Marquardt fit to the linear terms for each 
sinusoid frequency in a grid. For each frequency, a fit statistic was computed, given by 

A, .2 ,2 |0i - 0model(^i)P /n r\ 

i=l » 

where Xp^jy is computed for a fourth order polynomial model (no sinusoid) and 0j are the 
Crab phase residuals with errors Uj. 

First sinusoid frequencies in a wide range (0.5 to 4.0 cycles per 1000 days) were searched. 
This search found a peak in Ax^ at / 1.76 cycles/kday (equivalent to a period of 568 days) 
and a gradual increase in Ax^ as the frequency decreased below about 1 cycle/kday. Next 
a narrower search (1.26-2.26 cycles/kday), centered on the peak was performed to refine the 
flt. The best flt with frequency was / ^ 1.762 cycles/kday (or P fa 567.5 days) with a 
sinusoid amphtude {A^ + B'^y^^ pa 0.26 cycles. 

To estimate the probability of flnding a peak at this frequency by chance and to place 
errors on the parameters, two sets of Monte-Carlo simulations were performed. In the flrst 
set of simulations, the model above was flttcd to 10000 red-noise realizations and the bcst- 
flt frequencies were retained. The red noise consisted of simulated power-law noise using 
the time-domain method described in Appendix A with a power-law index m = 3 and 
noise strength 5" = 1.0 x 10~® cycles^/day^, equivalent to that found for the Crab pulsar. 
The same data sampling as the Jodrell Bank data was used. A distribution of the best-fit 
frequencies was obtained that had a maximum near the frequency 0.88 cycles/kday and a 
broad low frequency tail. In the second set of simulations a sinusoid with the parameters in 
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Table 2 was included. The resulting distribution of best-fit frequencies was tightly clustered 
around a maximum frequency of 1.76 cycles/kday. For the case with no sinusoid present, 
a total of 69 trials fell in the frequency range 1.736 to 1.790 cycles/kday, the full width of 
the distribution of best-fit frequencies with the sinusoid present, suggesting a probability of 
chance occurrence of 6.9 x 10~^. However, the largest sinusoid amplitude measured in this 
frequency range was ~ 0.11, about 2.3 times smaller than the amplitude measured in the 
actual Crab data. Also, the largest Ax^ = 6.9 x 10^, was about 5.4 times smaller than that 
measured in the data. These additional factors suggest that the probabihty that this peak 
is due to some systematic effect is smaller than 6.9 x 10~^. 

The errors on the fitted sinusoid parameters were estimated from the set of 10000 fits to 
the simulated red- noise realizations with a sinusoid present. The obtained parameter distri- 
butions were approximately gaussian and parameter errors were determined from gaussian 
fits to the distributions for the sinusoid frequency, cosine amplitude, and sine amplitude 
respectively. The parameters and their errors were also converted to their equivalent rep- 
resentation for a single cosine and phase^. Table 2 lists the best fit parameters along with 
their estimated errors. 

The low frequency oscillation in the 1982-1989 Crab phase residuals appears to be nearly 
periodic. This is suggested by the fact that when the power spectrum is recomputed after 
removal of the best fit sinusoid, essentially all excess power in the frequency range of the 
peak is removed. However, since less than five cycles arc present in the stretch of Jodrcll 
Bank data analysed here, and longer data stretches are interrupted by the occurrence of data 
gaps or large glitches, the long term coherence of the low frequency oscillation is not known. 



8.3. The 1986 Glitch 

We investigated the effect of the 1986 glitch (Lyne, Pritchard & Smith 1993) on the 
power spectrum. We tried simply deleting the glitch from the phase residuals or fitting a 
function to the glitch in order to smooth over the glitch interval. We excised the glitch over 
the interval TJD 6664.4 to TJD 6746.0. When this was done the spline interpolation filled in 
the empty gap with a linear segment. The computed window power spectrum (using a Hann 
window with a — 1) exhibited a slight change in slope from —3.09 ± 0.05 to —3.12 ± 0.06 
over the analysis frequency range 0.003 — 0.1 cycles/day. The mean power in flattened 
power spectrum dropped significantly from a level of 9.26(±0.49) x 10~^ to 5.3(±0.3) x 10~^ 
cycles^ /day^. The power at highest and lowest frequency portions of the PDS was unaffected. 



^0Qp(i ~ to) = ^0 cos(27r(t — to)/P — $o) where the single cosine amphtude is given by Aq = {A^ + B"^) i . 
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In contrast, when random sections of data of equal duration were removed, there was almost 
no effect on either the power spectral slope or the mean power level in the flattened spectrum. 

The second method used a fit to the glitch involving several steps. First the phase 
residuals from the preglitch interval one month prior to the glitch were fit with a linear 
function. This was extrapolated into the ghtch interval and subtracted. A glitch function of 
the form: 

(j){t) = 86400i/iTi(exp(-^/ri) - 1) 

where ui = 1.73 x 10~^ Hz and Ti = 6.3 days was then subtracted. A rising phase residual 
was left that was fitted with a quadratic polynomial and subtracted as well, leaving flat 
phase residuals with no systematic trends. When this set of functions was removed from the 
phase residuals, the phase residuals exhibited a linear trend in this interval. The result on 
the computed power spectrum was essentially the same as in the case of the simple excising 
above: a slope of —3.11 ± 0.06 and mean power level in the flattened power spectrum of 
5.6(±0.3) X 10~^ cycles^/day^ over the analysis frequency range 0.003 — 0.1 cycles/day. 

Either method of treating the glitch showed that a considerable amount of power resulted 
from the glitch but without any change in the power-law slope. Thus the noise strength 
estimated for the timing noise may be too high by a factor of two. 

9. Deeter Polynomial Power Spectral Estimation 

We computed power spectra of red PLN processes and of the Crab timing residuals 
using the Deeter polynomial method and compared them to power spectra computed us- 
ing windowing methods. We did this as a cross check on both power spectral estimation 
procedures and as a direct comparison with previously reported Crab timing analysis. The 
Deeter polynomial method of power spectrum estimation is directly applicable in the case of 
nonuniformly sampled data, unlike windowing methods, but produces power spectra of low 
spectral resolution. After comparing a number of procedures for estimating the power spec- 
trum of a red-noise process, Deeter (1984) concluded that a technique based upon sampling 
with orthonormal polynomials on a hierarchy of time scales produces a good power spectrum 
with modest frequency resolution. Specifically, Deeter (1984) proposed constructing power 
spectra with approximately octave frequency resolution by iteratively partitioning a data set 
into increasing numbers of shorter segments. For each segment, an amplitude is obtained by 
constructing a set of polynomials that are orthonormal on the segment and using the highest 
order polynomial from the set to "sample" the data by forming the sum of the product of 
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the time series with the highest order polynomial, i.e., 

a = ^(pkhk, (16) 

k 

where (fik is the measured phase and hk = p{tk) / cr^, where p{tk) is the (highest order orthonor- 
mal) polynomial evaluated at time tk, the time of the kth phase, and ak is the measurement 
uncertainty of the kth phase. The square of the amplitude, a, is an estimate of the power in 
the segment. For each partitioning, the average of the powers in the segments is an estimate 
of the power in the time series, at the time scale of the segments. By repeating the procedure 
for each partitioning, and by converting time scales to frequencies, a power spectrum of the 
white noise underlying the red noise is constructed. 

For a given set of times and weights, there is a unique set of orthonormal polynomials. 
For example, for the continuous interval [—1, 1] with uniform weight, they are the Legendre 
polynomials. It is not difficult to show that the amplitude obtained by sampling with the 
highest order polynomial in a set of orthonormal polynomials is equal to the coefficient of the 
highest order term in a polynomial fit to the data with the same order, /, as the highest order 
orthonormal polynomial. We chose to fit, rather than sample, and refer to this technique 
as the Deeter polynomial power spectrum estimation method or more simply as the Deeter 
polynomial method. We have applied this technique in the case of accretion-powered pulsars 
in Bildsten et al. 1997. 

As a direct comparison of the Deeter polynomial method to the windowing method we 
applied each method to the same simulated red PLN time-series. In figure 8 we show the 
average power spectra computed using each technique for the same 32 uniformly sampled 
noise reahzations. The noise is 1//^ PLN generated using the time domain method and 
detrended by a cubic polynomial fit. The spectra have been flattened for the case of a 
1//^ power-law, which involves dividing the Deeter polynomial power spectrum by a factor 
27r/. Sinusoids have also been added in four of the five spectra displayed. In general, 
when pure red PLN was present the two techniques gave consistent results as can be seen 
in the third spectra from the top, although the normalization of the Deeter polynomial 
method has to be adjusted somewhat depending on the assumed noise present. When strong 
discrete features were present in the power spectra, considerable power leakage occurred due 
to the low spectral resolution of the Deeter polynomial estimators and hence distorted the 
power spectra. In general, the presence of significant discrete features can cause considerable 
ambiguity in the interpretation of the Deeter polynomial power spectra. 

A Deeter polynomial power spectra was computed for the Jodrell Bank Crab timing 
residuals and is displayed in figure 1. The Crab timing residuals used in the computation 
are shown in the top panel of figure 1 and were subjected only to the dispersion measure 
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corrections and removal of outliers described in section 3. No averaging or interpolation 
was done. The computed power spectrum for the pulse frequency rate of the phases is 
shown in the middle panel of figure 1. The power contributed by measurement noise is 
shown by the triangles and has not been substractcd off from the total power shown by the 
solid circles, which is the convention used in this paper. If the timing residuals consisted of 
pure 1//^ noise in pulse phase as claimed by several previous analysis the power spectrum 
would be flat. A distinct rise with frequency can be observed at frequencies greater than log 
/ = —2.5, although the large amount of power at log / = —2.9 might cause one to argue 
that the spectrum is actually flat. When the spectrum is divided by 27r/ (bottom panel 
of flgure 1), which will produce a flat spectrum if 1//^ noise is present, a flat spectrum is 
indeed produced except for a peak near log / = —2.9. For comparison, we also plot the 
power spectrum (histogram) computed using a Hann window with a = 2 for the uniformly 
sampled time series displayed in the top panel of figure 6. The power spectra are similar 
except for a deficit in power in the Deeter polynomial spectrum near log / = —2.6, relative 
to the windowed spectrum. 

Since past analysis of the Princeton optical timing observations have concluded that the 
timing noise in the Crab pulse phase is 1//^ noise (Groth 1975c, Cordes 1980, Deeter 1981) 
rather than the 1//^ noise found above, we also investigated this data set to see why this 
was the case. In figure 9b we show a Deeter polynomial power spectrum calculated using 
the same Princeton optical data set used by Deeter (1981), which is displayed in flgure 9a. 
The power spectrum calculated by Deeter can be found as flgure 8 in Alpar, Nandkumar & 
Pines (1986). The spectrum is calculated for the pulse frequency derivative and hence will 
be fiat for pure 1//^ noise in pulse phase. The spectrum in our figure 9 is relatively flat 
but shows an increase in power at low frequencies. When flattened under the assumption of 
1 / p noise we see that the spectrum is again relatively fiat with a relatively greater excess of 
power at low frequencies very similar to that found for the JodrcII Bank data set. However, 
as the Princeton data set is somewhat shorter than the Jodrell Bank data set, the excess low 
frequency power appears more as a curvature in the spectrum rather than a discrete peak. 
We therefore conclude that the spectrum is consistent with either an interpretation as 1//^ 
noise in pulse phase or as 1//^ in pulse phase combined with a low frequency periodic or 
quasiperiodic process. In the absence of any motivation to assume that a large amplitude 
low frequency quasiperiodic process was present, clearly it would be more natural to assume 
that only 1//^ noise in pulse phase was present, although Groth (1975c) did speculate about 
the presence of a small contribution from 1//^ noise in pulse phase. 

A cubic spline interpolation of the Princeton optical timing residuals onto a uniform time 
series was made and a windowed power spectrum computed for the sake of completeness. The 
interpolated time series is shown in figure 9a offset below the actual timing measurements. 
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When a Hann windowed power spectrum with a — 2 was computed, a slope of —3.2 ± 
0.08 between the analysis frequency range of 0.003 — 0.1 cycles/day. The power spectrum 
flattened under the assumption of 1//^ noise is shown as the histogram in figure 9c. For the 
same frequency range, a mean noise level of (4.9 ± 0.37) x 10^^ cycles^/day^ is measured, 
about half that for the Jodrell Bank phase residuals. The windowed power spectrum is 
remarkably similar to the Deeter polynomial power spectrum and shows a similar large low 
frequency bump. The interpolated time series is clearly in error from the true time series 
since the gaps tend to be filled with linear segments. However, neither method of power 
spectrum computation has information within the gaps and the lowest order polynomials 
in the Deeter polynomial power spectrum also involve fitting across gaps and hence also 
involve a form of interpolation. Since the data are likely to be relatively smooth at the 
lowest analysis frequencies, these interpolations may not introduce any major distortions to 
the power spectrum. 



10. Discussion 

The timing noise in the Crab pulsar has been claimed to be either a white torque noise 
process or a quasiperiodic process. As a result of the analysis presented here, the timing 
noise or rotational irregularities in the pulse phase of the Crab pulsar appear to be composed 
of two components: a noise process with a 1//^ power spectrum and a periodic or nearly 
periodic process with a period of approximately 568 ± 10 days. A noise process with a 1//^ 
power spectrum in pulse phase becomes a process with a spectrum rising as / in the pulse 
frequency derivative, so the timing noise in the pulse frequency derivative (or torque for the 
case of a constant moment of inertia for the neutron star) has a blue power spectrum. 

The argument for 1//^ noise rests on the windowed power spectrum of the timing 
residuals obtained from the Jodrell Bank radio observations over the period 1982 - 1989. The 
windowed method of power spectrum computation was thoroughly checked with simulated 
noise data and appears to be a quite reliable method in this situation of relatively smooth 
residuals and high density sampling. The result is reinforced by a Deeter polynomial power 
spectrum of this same data which yielded a consistent result. A Deeter polynomial power 
spectral analysis of the Princeton Crab optical timing data previously analyzed by several 
investigators (Groth 1975b, c, Cordes 1980 and Deeter 1981), showed that this shorter and 
more sparsely sampled data set was consistent with either a pure 1/ noise process or a 
1/ P noise process with a large excess of low frequency power, thus explaining why previous 
analyses chose the simpler interpretation of pure 1//^ noise. 

Many theoretical studies have proposed models for the timing noise observed in radio 
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pulsars involving internal or external torques acting on the neutron star crust. A good 
summary and references can be found in D'Alessandro (1997). Many of these models have 
been developed to the point that they make predictions about the form of the power spectrum 
of the torque noise. In the vortex creep model (see Alpar, Nandkumar & Pines (1986)), 
several predicted torque power spectra are proposed that arc flat with a transition over 
some characteristic frequency range to another flat region or are flat at high frequencies but 
transition over to a 1//^ power-law at low frequencies. None of these models match the 
spectrum observed for the Crab pulsar. In the internal torque noise model of Jones (1990), 
the predicted power for the torque noise is a superposition of red and white components and 
also does not match the spectrum observed for the Crab pulsar. The external torque noise 
model of Cheng (1987) predicts either a blue or white spectrum for the torque noise but the 
blue component is a power-law rather than /. There does not seem to be any current 
model that predicts a blue pulse frequency derivative noise power spectrum rising as /. 

An excess of power in both the windowed and Deeter polynomial power spectrum is 
apparent at low frequencies indicating the additional presence of a strong quasiperiodic or 
periodic process. A known property of red power-law noise processes (or any noise process 
with a red spectrum) is that polynomial detrending will produce quasiperiodic residuals 
that are purely an artifact of the detrending. These quasipcriodicities have occasionally 
been considered real by investigators unfamiliar with this property of red noise. However, 
we show the quasiperiodicity present in Jodrell Bank Crab timing residuals is real. This is 
indicated by the presence of a peak in the low frequency power in the flattened windowed 
power spectrum rather than the deficit of power that would occur if the quasiperiodicity 
was spurious. The sharpness of the low frequency peak also argues against this being the 
result of a step in the power spectrum coupled with low frequency power removal caused by 
detrending although we did not investigate this possibility in depth. 

The coherence of the low frequency quasiperiodic oscillation was investigated by fitting a 
sinsuoid -|- quartic polynomial to the timing residuals using a linear least squares technique. 
The removal of this sinusoid also removed most of the excess power in the peak in the 
flattened windowed power spectrum (as well as in the Deeter polynomial power spectrum) 
indicating that the process is more nearly periodic than quasiperiodic. Another indication of 
the near periodicity is the comparison of the fitted sine to the 1st derivative of the Crab pulse 
phase timing residuals which revealed four clear oscillations after TJD 5500 (see figure 7). 
Before TJD 5400 a possible phase shift and amplitude change may have occured showing that 
the oscillation process may be incoherent on time scales longer than 2000 days. Inspection 
of the 1st derivative plot also revealed a possible small glitch occuring near TJD 5243 based 
on the occurence of a large amplitude spike that resembles the large amplitude spike caused 
by the known small glitch at TJD 6664. However, the data sampling around this second 
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possible small glitch is too sparse to make this conclusive (see figure 7). The occurrence of 
the glitch(s) does not seem to be related to the phase change in the long term oscillation 
since the phase was not clearly afi^ected by the ghtch at TJD 6664. 

The possible origin of the periodic or nearly periodic process has been discussed before 
by Lyne, Pritchard & Smith (1988). Explanations include oscillations of the superfiuid 
vortex lattice (Tkachenko oscillations), free precession of the neutron star or the presence of 
an orbiting companion. The 568 day period is roughly comparable to the 1000 day timescales 
observed in a few radio pulsars that may be exhibiting free precession (e.g Stairs, Lyne & 
Shermar, 2000). If the neutron star is freely precessing then long period pulse shape changes 
should be observable as well but no significant detections of this effect have been reported 
to date. If an orbiting companion is causing the observed rotation oscillations we can derive 
limits on orbital parameters. Assuming a binary companion with a circular orbit, we can 
determine an ax sin i — (8.65 ±1.0) x 10~^ It — sec from the cosine amplitude given in table 
2 and a mass function 



If we assume a Crab pulsar mass of = 1.4 Mq and a value of sini = 1 we can determine 
a lower limit on the companion mass of Mc > 3.2M0 and a projected orbital radius Oc sin i — 
750 ± 93 It-sec or 1.50 astronomical units. The presence of a planet around the Crab pulsar 
is an exciting possibility but requires confirmation by detection of the coherent oscillation on 
longer timescales before this could be considered the definitive explanation. Caution must 
be emphasized here as planets have often been claimed to explain the polynomial-like timing 
residuals of pulsars but have later failed to be confirmed (e.g. Konacki ct al. 1999), although 
it must also be emphasized that the longtcrm periodicity or quasiperiodicity here detected 
in the Crab pulsar is a real timing component and not merely an artifact of fitting red noise. 

Lastly, we speculate on a possible external origin for most of the power in the 1//^ 
timing noise component of the pulse phase. The mechanism causing the timing noise would 
have to have an amplitude that increases with analysis frequency, a characteristic shared 
by the gravitational perturbations caused by a body in a Keplerian orbit on its companion. 
Given that a planet may be in orbit around the Crab pulsar could there also be a disk 
of smaller bodies or planetesimals as well? The collective gravitational perturbations from 
a "disk" of planetesimals acting on the neutron star will be a source of timing noise. A 
planetesimal of mass m in a circular orbit with radius r around a neutron star with a mass 
Mj. — 1.4Mq will have a Keplerian orbital frequency (in cycles/day) given by: 



fiM) 



47r(aa;Sini)^ (Mc sin i)^ 



(6.8 ±2.4) X 10"^^ M, 



(17) 



Porb = 36.17r 2 
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where r is given in It — sec. The amplitude of the sinusoidal timing signature (in seconds) 
will be given by r = (m/Mj.)(r/c) where is the speed of light. Now if we assume that the 
average planetesimal mass scales with orbital radius, as might be the case if the planetesimals 
condensed out of a uniform density disk, then m — kr. Substituting into the relationship 
for r and solving for r in terms of Ugrb we get: 

T = 119.7(A;/M^)(i/<,,6)-t. (19) 

The superposition of timing signatures from many planetesimals over a wide range of orbits 
will result in a power spectrum of the pulse phase with a slope of —8/3 - very close to the 
observed slope of —3. 

We tested this idea using a Monte Carlo simulation by simulating a set of 128 plan- 
etesimals with random initial orbital phases and a range of orbital frequencies from 0.001 — 
1 cycles/day distributed uniformly in the logarithm. We fixed the constant k by choosing a 
planetesimal mass at r = 750 It — sec. For mass distributions with m — kr, the power spec- 
tral slopes tended to be near —3.7, too small to fit the Crab pulsar. For a planetesimal mass 
distribution with m — kr^ and a planetesimal mass of 1.62M0 at r = 750 It — sec, the average 
slope of the power spectrum was near —3 over the frequency range 0.003—0.1 cycles/day, with 
an average noise level in the flattened power spectrum of 1 x 10~^ cycles^ /day^ - matching 
that observed for the Crab pulsar. To fully model the Crab timing residuals the planetesi- 
mal at r = 750 It — sec must be replaced with a planet of mass M > 3.2M0. Whether such 
a planetesimal mass distribution and set of orbits is reasonable and physical is a question 
beyond the scope of this paper. 

The power spectrum of the timing caused by a planetesimal disk can be calculated 
analytically as well. The phase signature of the the disk of planetesimals is 

A0(i) = 5] cos(20i/fet + ^k) (20) 

k 

where the sum is over plantesimals k, and Xk, i^k, and ipk are the phase amplitude, orbital 
frequency, and orbital phase respectively of plantesimal k. The phase amplitude Xk is given 

by 

^k — J^oT^k sin i mkM~^c~^ (21) 

with Tfe the orbital radius and the mass of plantesimal k, i is the inclination of the disk, 
and M the total system mass. The orbital radius is related to the orbital frequency by 



r{u) = (G'M)^/3^27rz/)-3 



(22) 
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If we assume the orbital phases xjjk are independent then the autocorrelation function of A(f) 
is 

^A0(t) ^k<J2^k cos{2nUkT) > . (23) 

k 

If the density of plantesimals per unit radius r is given by D{r) and the plantesimal mass 
distribution at radius r is given by P{m\r) then the autocorrelation function is expected to 
be 

S{iy) cos(27ri/r)(ii/ (24) 

where we have changed integration variables using equation 22, and 

S{u) = Inu^ sin^ iGM-^c-^{27ru)-^D{r(u))m^(r{u)) (25) 

with 



m (r) = J P{m\r)dm (26) 

being the mean square plantesimal mass are radius r. 

Slu) is the power spectrum of the phase component A0. If D{r)m^{r) is a powerlaw in 
r with index a, then <S(i/) will be a powerlaw in u with index — (3 + 2q;/3). The measured 
power spectra implies that a ~ 0. Prom the strength of the power spectrum we find 

that 

D{r)m^{r) w 2.1 sin'^ iM^Au"^ (27) 
in the radius range 0.13 < r < 1.3Au. 

Remnant disks of heavy elements formed by fallback material from the supernova explo- 
sion that gave birth to the neutron star have been proposed as explanations for phenonmena 
observed in radio pulsars and anomalous X-ray pulsars (e.g. Lin et al. 1991, Alpar 2001). A 
potential planet or planets around the Crab pulsar would probably have to come from preex- 
isting planets that survived the supernova rather than being assembled from planetesimals 
given the short time (< 1000 years) since the Crab supernova occurred although the effect 
of a pulsar wind on a dust disk might significantly speed up the formation of planetesimals 
by concentrating the orbiting dust into preplanetesimal pockets. In addition to the cases 
of a pure fallback disk or surviving planets one can also envision mixed scenarios in which 
supernova surviving planets interact with a supernova fallback disk or a planetesimal disk 
forms through the collision of a planet with another body. 
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The presence of a planetesimal disk causing the timing noise in the Crab pulsar and 
possibly in other radio pulsars with significant timing noise would have some wider implica- 
tions if validated. The fall back scenarios proposed for anomalous X-ray pulsars (AXP's, for 
a review see Mereghetti et al. 2002) have tended to consider only gaseous accretion disks. 
Large gaseous disks could supply the spin-down torque observed in these systems and supply 
accreting material to produce the persistant X-ray emission, but cannot explain the recently 
observed X-ray bursts (Gavriil et al. 2002) that support the idea that the related X-ray 
burst emitting soft gamma repeaters and anomalous X-ray pulsars form a single class of 
objects. In addition, searches for optical counterparts do not support the presence of large 
gaseous disks in such systems (e.g. HuUeman et al. 2000). Thus, by default, the alternative 
magnetar theory has gained much recent support. In contrast to a large gaseous disk, a 
toroidal swarm of planetcsimals would be more easily hidden from view and could supply 
material in the form of gas for steady accretion as planetcsimals are tidally broken up near 
the neutron star and form an inner gaseous ring or as dust formed further out, perhaps 
through plantesimal intradisk collisions, migrates inwards under the Poynting-Robertson ef- 
fect or from drag in a pulsar wind. Planetesimals beyond the tidal radius could also supply 
material for series of occasional X-ray bursts if a planetesimal gets fragmented as a result of 
an intradisk collision and the low angular momentum pieces gradually rain down onto the 
neutron star in an event analogous to a meteor storm. Alternatively, a planetesimal could 
have its orbit disrupted into a comet-like path by a close encounter with one of the largest 
disk members and eventually collide with the neutron star causing a giant outburst. Series 
of smaller outbursts could also occur before, during and after the main body collides with 
the neutron star, if the plantesimal is broken up by a close orbital pass by the neutron star 
prior to or during the orbital pass of the primary coUision. Either scenario would cause 
occasional bursts of X-rays to be emitted by the neutron star. Thus a plantesimal disk has 
the potential of explaining both the persistant X-ray emission observed, the slowing down of 
the neutron star spin, sets of occasional small and giant X-ray bursts and the observed lack 
of a large gaseous accretion disk. In addition, the observations of transient infrared emission 
associated with several AXP's after X-ray bursts have been detected (e.g. Israel et al. 2002) 
could be related to the formation of these temporary orbiting dust or debris clouds. 

Chandra X-ray observations reveal an equatorial torus of material around the Crab 
pulsar (Weisskopf et al. 2000) that might be supplied by ablation and/or tidal disruption 
of plantesimals from an inner disk. The radius of the observed torus is on the order of 0.4 
parsecs, while the disk considered here is within 1.5 astronomical units of the pulsar - far 
too small to be imaged directly. 

The fast rotation of the Crab pulsar would place it in a propellor phase in which the 
outward flow of momentum from the pulsar exceeds the gravitational energy of the incoming 
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material. This would result in a particle outflow and possibly cause the observed jet from 
the Crab pulsar. The rate of spindown torque on the neutron star would also be increased 
relative to the vacuum case, lowering the braking index from the canonical value of 3 for 
the emission of pure magnetic dipole radiation. However, to observe gravitational effects on 
the neutron star it is not necessary that the planetesimals or any derivative dust must be 
interacting directly with the pulsar magnetosphere at the present time. 

A group of planetesimals with a total mass of a few dozen Earth masses could plausibly 
cause the timing noise in the Crab pulsar and would have a mass similar to that required of 
the postulated "fallback" disks of the anomalous X-ray pulsars and soft gamma repeaters. 
If the scenario outlined here is correct, the Crab pulsar could eventually evolve into an 
anomalous X-ray pulsar or soft gamma repeater when the neutron star rotation slows to the 
point that accretion becomes possible. If other pulsars have planetesimal disks, then the 
power spectral indices of their timing noise might be expected to have a fairly wide range, 
caused by differing disk mass distributions and one might observe real quasiperiodicities as 
well in their timing noise due to the largest bodies orbiting the neutron star. 

We thank Andrew Lyne for kindly providing the Jodrell Bank Crab timing data and 
Robert Pritchard for help in the initial analysis. We also thank Paul Boynton for useful 
discussions and for help in the early stages of this project. 

A. Generation of Red Power-Law Noise Using Power-Law Weighting of White 



Power-law noise (PLN) processes with negative even integer spectral indices can be 
produced by the repeated integration of a white noise process e{t). For example, a random 
walk, with a 1//^ power spectrum, is defined as: 



We consider e{t) to be a zero-mean white noise process defined for all t with an autocovariance 
function < e(t)e(t + r) >= a'^S{T), where S{t) is the Dirac delta-function. An equivalent way 
of defining r2{t) is as the convolution of white noise with the unit step function, H{t): 



where H{t) is the standard Heaviside step function defined as: H{t) — 1; t > and 
H{t) — 0; t < . Now a 1//^ noise process, produced by the integration of r2{t), can be 
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made instead by convolving H{t) with r2{t), but this is equivalent to first convolving H{t) 
with itself and then convolving the resulting function with e{t): 

H{t) * {H{t) * e{t)) = {H{t) * H{t)) * e{t). (A3) 

The self convolution of H{t) results in a ramp function: tH{t). Therefore r4{t) can be 
generated directly as the convolution of tH{t) with white noise: 

r^it) = {tH{t)) * e{t) = j {t- t)e{t)dt. (A4) 

J —oo 

The convolution of e{t) with t^H{t) will generate 1//^ noise. In general, even-integer red 
PLN can be generated directly from white noise as: 

rm{t) = ^ {t- tfe{t)dt (A5) 

where the power-law spectral index m = 2(A; -|- 1) and m — 2, 4, 6 . . . without the need for 
repeated integration. Note that a normalization factor of ^ has been added. Given this 
formulation for r„j(t), there is no a priori reason to restrict k to positive integer values and 
we consider the generalization of this formulation to noninteger k. 

The noise process described in equation (A5) is a member of the class of general linear 
noise processes (see e.g. Priestley, 1981, p. 177) with the form: 

r{t) = / w{t- t)e{t)dt (A6) 

J —oo 

in which a weight function of the form w{t)H{t) is convolved with white noise. This form 
has many easily analyzable properties. For example, the ensemble averaged PDS can be 
calculated from the multiplication of the Fourier transform of w{t)H{t) with that of e(t). 
Thus r{t) can be described in terms of the "noise strength" S and the weight function of the 
form w{t)H{t) in the time domain or the spectral shape of the PDS in the frequency domain. 
In the case of discrete white noise, the weight function can be viewed simply as a "step" 
so that the noise process can be viewed as an accumulation of random steps with a similar 
form. For example, a random walk is an accumulation of the step functions H(t), occuring 
with random amplitudes and at random times. If the duration of the step is finite the noise 
process will be asymptotically stationary since two values of the noise process separated by 
a sufficiently large time interval will be uncorrelated. In the well known case of "shot noise" 
the step is a decaying exponential, so values of the noise process are correlated on short 
timescales and approach an uncorrelated state on long timescales. This is reflected in the 
PDS of shot noise, which is fiat at low frequencies but turns over to a power-law at high 
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frequencies. For red power-law noise proceses a correlation exists between a given value and 
all past values of the time-series, but the statistical properties are the same at any time t so 
the process is stationary. 

Let us look at a specific power-law weight function before moving onto the more general 
case, lik — so w{t)H{t) — -^t~^^'^H{t), where the generalized definition of k\ — r(l-|-A;) 
has been used, then the Fourier transform will be given by: 

1 f~^°° 

F{f)^— t-^/^{cos{2nft)-ism{2nft))dt. (A7) 
VTT Jo 

The resulting Fourier transform is: 

n/) = (l-0(:r^) (A8) 
with a power density spectrum of the form: 

P{f) = F{f)F*{f) = ^. (A9) 

Multiplication by the PDS of e{t) creates the well known 1// noise power spectrum. 

The general power-law problem involves solving Fourier transforms of the weight func- 
tion W{t) = '^H{ty. 

F{f) = T7 / t^{cos{27r ft) - ism{2Tr ft))dt. (AlO) 
^' Jo 

The Fourier integral will not in general converge since the time-series has infinite "en- 
ergy" but we can solve directly for a few special cases. For example, if — 1 < /c < the 
integrals have the solutions: 

and the power spectrum is given by: 

^^^^ ^ Wmm + W (sin(M^)2cos(Mzi)2) (2^/)"'^'^'^- (A12) 

When k — then (All) and (A12) reduce to the previous results (A8) and (A9) for 1/ f 
noise. 
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Although the Fourier transform can't be solved directly it can be solved as a limit. If 
we add an exponential cutoff we get an integral of the form: 

1 r+oo 

F{f)^— e-U''{cos{2nft)-ism{2nft))dt (A13) 
^' Jo 

The integral solutions to the cosine and sine portions of the integral can be found in integral 
tables and are, respectively: 



1 r 



2(j, + (27r/)2)'+i ^ ' 

and 

In the limit as A — > 00 the Fourier transform will approach the value: 

= (-^)^'+i(27r/)-('^+^) (A16) 
and the power density spectrum will approach the power-law form: 

P(/) = (27r/)-^(^+^). (A17) 
The full ensemble averaged PDS will be given by: 

2 

P(/) = (A18) 

where a1 is the mean power level of the underlying white noise (equal to the noise strength 
S) and m = 2(A; + 1) is the power-law spectral index. Now equation A5 can be viewed as 
a general formula for producing power- law noise by the convolution of a step t^Hif) with a 
white noise process. 



A.l. Time Domain Properties 

In practice, any physically realizable red noise time series will first be observed at a time 
t = and have started at this or some earlier time. The linearity of the convolution integral 
allows the time series to be divided into the form: 

r™(t) = -^ / {t-tYe{t)dt + ^ f {t-t fe{t)dt (A19) 



-34- 



where T can range to — oo and m — 2{k The first integral referring to the past will 

contribute a smooth function to the time series after t = 0. Note that in this form, Tmit) is 
divided into two nonstationary processes. When e{t) is viewed as a series of discrete impulses 
the time series rm{t) can be viewed as an accumulation of "steps" of the form s{t) = t^H{t). 
The smooth function will have a form related to that of an individual step. If the index k is 
a positive integer, the smooth function will have the form of a simple polynomial. For the 
random walk, k=0 (1//^ noise), the time series consists of a series of randomly ascending or 
descending "stair steps" of variable amplitude and duration. The smooth function will be a 
constant: 



e{t)dt = Co. (A20) 



IT 



For k=l (1// noise), the time series consists of a series of connected linear segments 
( "ramps" ) and the smooth function will be a linear trend: 



' [t - t)e{t)dt = t / e{t)dt - I te{t)dt = Cot - Ci 
T Jt Jt 



(A21) 



and for k — 2 (1//^ noise), the time series consists of quadratically shaped segments and 
the smooth function will be a quadratic trend and so on. Note that a polynomial fit of an 
order I > k will remove the smooth function component completely when k is an integer. 

The ensemble averaged expectation and variance of power-law noise at a time t can be 
easily calculated. Each point of a power-law noise realization simply consists of a linear com- 
bination of white noise, so the statistics will be directly related to the white noise statistics. 
The expectation, < rm{t) > will be given by: 

< r^t) >= ^ (^f it - t)'e{t)dt'^ = {e{t)) J.fit- tfdt. (A22) 

If the times series is divided into two finite portions as in A19 and k ^ —1 the expectation 
has the simple form: 

<U^)>=^^(^ + 7^r^ (A23) 

A zero mean white noise process implies that < rjn{t) >= as well. 

The variance of r^(i) as a function of time, < rjn{ty > — < r^(t) >^ can be calculated 
from: 




J' {t-t)''e{t)dt^ 
f {t - tf{t - tf Ut)e{t)) dt'dt" 

D J —OO 
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^yj^jt-t'rdt (A24) 

where use is made of the fact that for white noise: 

(^e{t)e{t")'^ = a^J{t -t") (A25) 

where S{t) is the Dirac delta-function. Note that the variance < rm{tY > at time t is in fact 
infinite for all values of k. If the time series is divided into two finite portions as in A19 with 
zero mean white noise being used and k ^ —0.5, the variance is finite and reduces to the 
simple expression: 

Therefore, if one could measure < rm{ty > in this way, by averaging over a number of noise 
realizations for example, the PDS spectral index m = 2{k + 1), and the underlying white 
noise variance cr^ could be estimated in a straight forward manner. 

The autocovariance can also be derived in a similar manner. The autocovariance is, for 
< rm{t) >= 0, given by: 

<r^(iK(t + T)>= (l^y^l^f {t-tfe{t)dt^ (^j'^\t + r-t)'e{t)dt^y (A27) 

This expression reduces to: 

< rm{t)rrr^{t + t) >= (^)' (^j^^ (t - t)' (t + T - t )'dt^ (A28) 

with r > 0. Note that the autocovariance function depends on both the current time t as 
well as the lag r. The Fourier transform of the autocovariance function has the form: 

/oo 
e-''"^" < rm{t)rm{t + r) > dr (A29) 
-oo 

The exponential can be broken into the two terms e~^'^'^^'^ — e""*^'^-^^'^""^* -*))e-*27r/(t -t) 
after rearranging, the integral divides into two terms, one of the form AlO and one of its 
conjugate that can be solved in the limit as in A13. The Fourier transform will, in the limit, 
approach the power-law form given by A18. 

Expression A27, after removing the past smooth function, has the form: 

< rm{t)rm{t + r) (^)' (^j\t - tf{t r - tfdt^ (A30) 
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which reduces to A24 when r — 0. For the three cases in which k=0, 1 and 2, i.e. for 1//^, 
1//^ and 1//^ PLN this expression has the following forms: 
For k=0: 

< r2{t)r2{t + t) >^ a^^t (A31) 

for k=l: ^ 

< r4{t)r4{t + r) >= -a^t"^ {2t + 3t) (A32) 

and for k=2: 

< re(t)n(t + r) >= — (j'^' (6^' + ISir + lOr') (A33) 

Note the time dependence of the autocovariance for r4(t) and rQ{t) after removal of the past 
smooth function. When the substitutions t — and t + r — t> are used then A31, A32 
and A33 are seen to be the continuous white noise version of equations (23), (38) and (43) 
derived by Groth (1975b). 



A. 2. Discrete Red Power-law Noise 

Any realizable red PLN process, in addition to being of limited duration, will be com- 
posed of steps of a finite size occurring at a finite rate. We now consider a discrete model 
of red PLN that is useful for generation of actual samples with Monte-Carlo routines. A 
discrete form of (A5) can be obtained by replacing e{t) with the discrete white noise form: 

oo 

e{t) = Yl ^^■'^(^^•) (A34) 

j=-oo 

which produces the discrete power-law noise process: 

-| Jmax 

^mit) = ^ ^(t-t,)'i/(t-t,)6(i,) (A35) 

' j=-oo 

where the power-law spectral index m = 2{k + 1) and H{t) is the unit step function as 
defined previously. The last white noise "event" in the time-series, e(tj), occurs at a time 
^jmax such that tj^^^ < t < We choose the duration between events At = tj — tj-i 

to be a random variable with a probability density function: p{At) = Rge~^^^^, i.e. the 
"waiting time" distribution of the Poisson process. The t{tj) are drawn from a discrete zero 
mean gaussian white noise process with a variance 0"^ (At^) = Atpcr^ and a mean rate of step 
generation Rg = such that Atg —< At > (a noise process sometimes referred to as white 
Poisson impulse noise). The noise strength is given by S — Rga^{Atg) and along with the 
index k and Rg completely specify the type of noise process. Processes with the same noise 
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strength can be made by the frequent accumulation of many small steps or with fewer, larger 
steps in the same time interval. 

The physical units of S depend on the units of fm{t) and the specific power- law. For 
example, let f^(t) be noise in the rotation phase of a pulsar so the physical units are cycles 
and let time be given in seconds. Then the units of €{tj) will be cycles sec"''. Because 
a^(Atg) — the units of S are cycles^ sec"^^"'"^''^. If frn(t) is noise in the rotation 

frequency then the units of €{tj) will be Hz sec"'' and the units of S will be Hz^ sec"^-'"'"^''^. 
If noise in the rotation frequency derivative then S will have units of (Hz/sec)^ sec"^^''"^''^, 
etc. 

Any physical or simulated time-series will begin at a time that we can arbitrarily set at 
t — 0. Then f^(i) will have the form: 

-| Jmax 

^m(t) = ^ - ^^^'^(^ - *Mtj)- (A36) 
■ j=o 

In this form it can be seen that "Tmit) is simply a generalization of the discrete PLN models 
contained in Groth (1975b), in which k=0, l,and 2 correspond to his phase, frequency and 
slow-down noise definitions respectively. Wc can also rewrite e{tj) as a scaled version of a unit 
variance time-series e„(tj) such that: e{tj) = ^JSAtg€u{tj) where a'^^{Atg) = 1. Equation 
(A35) then becomes: 



/ Q A J- Jrnax 

it) = - *^)'^(^ - ^^■)^«(^^)- (A37) 

j=Q 



and is useful for generating discrete noise realizations of noise strength S. Now if f^(i) above 
is used in its entirety to represent a power-law noise realization, the low frequency power 
contributed by the past will be too small so for an accurate representation, a sufficiently 
long portion of the early time series needs to be discarded. 

Following Groth (1975b), we now derive the mean < fm(t) > and autocovariance < 
^m(^)^m(^ + t) > for equation (A36). Now < fm{t) > is given by: 



I Jmax \ 
-j^ Jmax 

= yJ2<it-tj)'H{t-tj)><e{tj)> 



j=0 

Jmax 
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f\t - tfdt 
Jo 



_ < e{t) > 

where < ^{Atg) > is the mean on a time-scale Atg = Rg^- This is the same equation as 
(A23) and when k=0, 1, and 2 reduces to Groth's results for phase, frequency and slow-down 
noise. 

The autocovariance is given by: 

+1 

< rm{t)Tm{t + r) > = ^ - t^fH{t - i,)e(i,) j E + ^ " *^)'^(* + ^ " ^^X^^) 

2 Jmax Jmax +1 



(-1 \ Z Jmax Jmaxii' 
m) {it-tjnt + r-t.rH{t-t,)H{t + T-t.)){e{tMtf)) 

j=0 j'=o 

V ■ / j=o 



where t^^^^ < t < tj,na.+i and t^Vn^.+z < t + r < tjVna.+/+i- The autocovariances derived 
by Groth (1975b) for his phase, frequency and slow-down noise definitions (k=0, 1 and 2 
respectively) can be obtained from equation (A39) by the substitutions i> = i -|- r in the 
second term in the integral and — t elsewhere. 



B. Generation of Red Power-Law Noise Using a Modified Timmer-Konig 

Method 

A method for simulating power-law noise realizations has been proposed by Timmer 
& Konig (1995). Real and imaginary values of the Fourier amphtudes corresponding to a 
set of discrete Fourier frequencies are generated from a zero mean gaussian distribution and 
then weighted by a chosen Fourier frequency dependent power-law. The amplitudes are then 
Fourier transformed into the time-domain to generate a uniformly sampled time-series. A 
normalization for a particular noise strength is not given in their paper. We present below, 
in a recipe form, a modification of this method that can be used to generate power-law noise 
realizations with noise strength 5*. 

To generate a power-law noise realization with N points we created two zero mean unit 
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variance gaussian white noise vectors Xj and yj, j — 0,1 . . . N — 1 and a corresponding 
weighting vector wj = (j)^ where m is the desired power-law index. For the conventions 
used below, we assume that N is an even integer. We then formed a complex vector of Fourier 
amplitudes r = A [0, (tUjXj, (w7v-i-ja:Ar-i_j, W7v-j-i|/7V-i-i)*] where the first member 

of a complex pair (wjXj, WjUj) is real and the second imaginary and the asterisk indicates a 
complex conjugate. The third set of terms contains the negative Fourier frequencies in which 
j = 1 . . . — 1, while for the second set of positive Fourier frequencies, j = . . . — 1, so 
that the complex vector r contains a total of 2A'" terms. The Nyquist term was made purely 
real so that Vn-i — 0. The normahzation constant A — {SN)^{^)^ and was chosen so that 
the generated realization has an average noise strength S. A realization was then generated 
as a time series via the discrete inverse Fourier transform: 

^ 2N-1 

^(^'^) = 27V ^ '■^''"^ A; = 0,l,...2Ar-l. (Bl) 

1=0 

While R{tk) is a time-series of N independent points, only the first N/2 should be used. 
The autocorrelation function of the times-series produced by equation Bl is approximately 
correct for lags less than N/2, but is wrong for larger lags. This is due to the cyclic nature 
of the discrete Fourier transform, which treats the last point as the nearest neighbor of the 
first. Use of the full time series results in an under representation of low frequency power. 
Further discussion can be found in section 7 of this paper. 
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Fig. 1. — Top: Crab arrival time residuals obtained from radio observations by Jodrell 
Bank after removal of a quartic polynomial and correcting for dispersion effects. Middle: 
Power spectrum of the frequency derivative computed from the above residuals using the 
Deeter polynomial method (solid circles). Measurement noise power is shown by the trian- 
gles. Bottom: Deeter polynomial power spectrum (solid circles) after dividing by 27r/. The 
measurement noise power (triangles) has been similarly corrected. The histogram displays 
the windowed power spectrum of the arrival time residuals displayed in the top panel using 
a Hann window with a — 1. 

Fig. 2. — Examples of discrete rcd-noisc time scries showing resolved steps. The number on 
the far right gives the power-law index of the process. The bottom time scries shows the 
discrete white noise realization from which the red-noise time series were generated. 

Fig. 3. — Sample of red noise realizations of increasing redness generated from the same 
underlying white noise realization (left panel set). The power-law index, m, for noise 
is given in the inset on each left panel set. The right panel set shows the same realizations 
after detrending by a cubic polynomial fit. 

Fig. 4. — Top: Example of a cubic polynomial detrended 1//^ noise reahzation with a noise 
strength S* = 1 x 10~^ cycles^/day^ (approximately equal to the Crab pulsar) generated 
using the time-domain method. Middle: Hann windowed [a = 2) power density spectrum 
(PDS) showing power-law slope of —2.998 ±0.077 from fit to power between vertical dashed 
lines. For comparison, the smoother curve offset above is an average PDS of 1024 identically 
generated realizations. Bottom: Flattened PDS after multiplication by (27r/)^, where f is 
frequency. The smoother average flattened PDS is also shown offset above where the highest 
frequency power exhibits an increase by 2x due to aliasing. 

Fig. 5. — Averages of 1024 flattened power spectra for power-law noise realizations of unit 
strength. A high frequency sinusoid (visible as peak near / = 0.1 and white noise have 
been added to each realization. Power-law spectral index m shown to left of the flat lines 
that display the predicted pure power-law form for unit noise strength. Power spectra with 
m > 2.0 have been offset in power for display purposes, otherwise scale is the same. The 
average spectra of the white noise alone (for / > 0.08) is shown as the sloping hues. The 
lower spectra of each pair show average power spectra after realizations have been detrended 
by cubic polynomials. 
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Table 1. Crab pulse phase ephemeris 



Parameter 



Value 



U{MJD) 
(f)o{cycles) 

MHz) 

z>o(10-i° Hz s- 
i>o(10-2° Hz s- 
uo (10-30 Hz s 



45395.94232638889 
4.445(6) 
30.04669082632(13) 
-3.80599313(13) 
1.214636(13) 
-1.425(16) 



Table 2. Fitted Cosine parameters 



Parameter 


Value 


to(MJD) 


46390.0 


A{cycles) 


0.23 ±0.03 


B{cydes) 


0.12 ±0.03 




1.09 ±0.116 


P{days) 


568 ± 10 


AQ{cydes) 


0.26 ±0.03 
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Fig. 6. — Top: Crab arrival time residuals obtained from radio observations by Jodrell Bank 

after removal of a quartic polynomial, correcting for dispersion effects and interpolating onto 
a uniform grid. Middle: Hann windowed {a = 1) power density spectrum showing power-law 
slope of — 3.09±0.05 from fit to power between vertical dashed lines (frequency range 0.003 — 
0.1 cycles/day). Horizontal line shows calculated power level of statistical measurement 
noise. Bottom: flattened PDS after multiplication by (27r/)^, where f is frequency. Note the 
highly significant peak at f — 0.0018 cycles/day. Mean power level (between vertical dashed 
fines) is 9.26(0.50) x 10~^ cycles^/day^. Sofid vertical line shows cosine frequency from table 
2. 

Fig. 7. — Top: The numerical 1st derivative of the uniformly binned spline fit to the Crab 
timing residuals. The superposed sinusoid is the derivative of the sinusoid fit to the Crab 
timing phase residuals. Middle: A closeup of the Crab timing residuals at the time of the 
small glitch occuring at JD - 2440000.5 = 6664. The raw mesurements are plotted in top, 
the spline interpolation underneath. Errors are plotted for the raw measurements but are 
generally smaller than the symbol size. Bottom: A similar closeup of the region associated 
with the first large excursion in the numerical 1st derivative. A possible glitch may have 
occured here. 

Fig. 8. — Comparison of average power spectra computed using the Deeter polynomial 
method (diamonds with error bars) and windowed power spectra using a Hann window with 
q; = 2 of the same 32 uniformly sampled noise realizations. The noise is 1//^ PLN generated 
using the time domain method and detrended by a cubic polynomial fit. The spectra have 
been fiattened for the case of a 1//^ power-law. High (top two spectra) and low (bottom 
three spectra) frequency sinusoids have been added to the noise realizations to show the 
effect of discrete features on the power spectrum estimation. The power spectra have been 
shifted in power for display purposes but the noise strength is the same in all cases. 

Fig. 9. — Top: Crab pulsar optical timing pulse phase residuals from Groth (1975c) for period 
after 1969 glitch. Offset below is a cubic spline interpolation onto a uniform time series. 
Middle: Solid circles mark the Deeter polynomial power spectrum of the pulse frequency 
derivative (ffat for 1//*^ noise in pulse phase). Triangles connected by dotted line mark 
computed measurement noise spectrum. Dashed line shows power level measured by Deeter 
(1981) in his analysis of these data. Bottom: Above PDS after multiplication by (27r/), 
which assumes that the real spectrum is 1//^ noise in pulse phase. Note the "bump" in the 
power near log frequency —3.0. Dashed horizontal line shows noise power level measured 
from Jodrell Bank measurements. Solid vertical line marks a period of 568 days. Histogram 
shows Hann windowed power {a — 2) power of spline interpolation flattened assuming an 
power-law index of —3. 
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Fig 1. Top: Crab arrival time residuals obtained from radio observations by 
Jodrell Bank after removal of a quartic polynomial and correcting for disper- 
sion effects. Middle: Power spectrum of the frequency derivative computed 
from the above residuals using the Deeter polynomial method (solid circles). 
Measurement noise power is shown by the triangles. Bottom: Deeter polyno- 
mial power spectrum (solid circles) after dividing by 27r/. The measurement 
noise power (triangles) has been similarly corrected. The histogram displays 
the windowed power spectrum of the arrival time residuals displayed in the 
top panel using a Hann window with a — 2. 
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Fig 2. Examples of discrete red-noise time series showing resolved steps. The 
number on the far right gives the power-law index of the process. The bottom 
time series shows the discrete white noise realization from which the red-noise 
time series were generated. 
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Fig 3. Sample of red noise realizations of increasing redness generated from 
the same underlying white noise realization (left panel set). The power-law 
index, m, for 1//"* noise is given in the inset on each left panel set. The right 
panel set shows the same realizations after detrending by a cubic polynomial 
fit. 
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Fig 4. Top: Example of a cubic polynomial detrended 1//^ noise realization 
with a noise strength 5" = 1 x 10~^ cycles^/day^ (approximately equal to 
the Crab pulsar) generated using the time-domain method. Middle: Hann 
windowed {a — 2) power density spectrum (PDS) showing power-law slope of 
— 2.998±0.077 from fit to power between vertical dashed lines. For comparison, 
the smoother curve offset above is an average PDS of 1024 identically generated 
realizations. Bottom: Flattened PDS after multiplication by (27r/)^, where f 
is frequency. The smoother average flattened PDS is also shown offset above 
where the highest frequency power exhibits an increase by 2 x due to aliasing. 
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Fig 5. Averages of 1024 flattened power spectra for power-law noise realiza- 
tions of unit strength. A high frequency sinusoid (visible as peak near / = 0.1 
and white noise have been added to each realization. Power-law spectral in- 
dex m shown to left of the flat lines that display the predicted pure power-law 
form for unit noise strength. Power spectra with m > 2.0 have been offset in 
power for display purposes, otherwise scale is the same. The average spectra 
of the white noise alone (for / > 0.08) is shown as the sloping lines. The lower 
spectra of each pair show average power spectra after realizations have been 
detrended by cubic polynomials. 
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Fig 6. Top: Crab arrival time residuals obtained from radio observations by 
Jodrell Bank after removal of a quartic polynomial, correcting for dispersion 
effects ano? interpolating onto a uniform grid. Middle: Hann windowed {a — 1) 
power density spectrum showing power-law slope of —3.09 ± 0.05 from fit to 
power between vertical dashed lines (frequency range 0.003 — 0.1 cycles/day). 
Horizontal line shows calculated power level of statistical measurement noise. 
Bottom: flattened PDS after multiplication by (27r/)^, where f is frequency. 
Note the highly significant peak at / = 0.0018 cycles/day. Mean power level 
(between vertical dashed hues) is 9.26(0.50) x 10"'' cycles^/day^. Solid vertical 
line shows cosine frequency from table 2. 
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Fig 7. Top: The numerical 1st derivative of the uniformly binned spline fit 
to the Crab timing residuals. The superposed sinusoid is the derivative of the 
sinusoid fit to the Crab timing phase residuals. Middle: A closeup of the Crab 
timing residuals at the time of the small ghtch occuring at JD - 2440000.5 = 
6664. The raw mesurements are plotted in top, underneath (diamonds) is the 
spline interpolation. Bottom: A similar closeup of the region associated with 
the first large excursion in the numerical 1st derivative. A possible glitch may 
have occured here. 
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Fig 8. Comparison of average power spectra computed using the Deeter poly- 
nomial method (diamonds with error bars) and windowed power spectra using 
a Hann window with a — 2 oi the same 32 uniformly sampled noise real- 
izations. The noise is 1//^ PLN generated using the time domain method 
and detrended by a cubic polynomial fit. The spectra have been flattened for 
the case of a 1//^ power-law. High (top two spectra) and low (bottom three 
spectra) frequency sinusoids have been added to the noise realizations to show 
the effect of discrete features on the power spectrum estimation. The power 
spectra have been shifted in power for display purposes but the noise strength 
is the same in all cases. 
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Fig 9. Top: Crab pulsar optical timing pulse phase residuals from Groth 
(1975c) for period after 1969 glitch. Offset below is a cubic spline interpolation 
onto a uniform time series. Middle: Solid circles mark the Deeter polynomial 
power spectrum of the pulse frequency derivative (fiat for 1/ noise in pulse 
phase) . Triangles connected by dotted line mark computed measurement noise 
spectrum. Dashed hne shows power level measured by Deeter (1981) in his 
analysis of these data. Bottom: Above PDS after multiplication by (27r/), 
which assumes that the real spectrum is 1/ p noise in pulse phase. Note the 
"bump" in the power near log frequency —3.0. Dashed horizontal line shows 
noise power level measured from Jodrell Bank measurements. Solid vertical 
line marks a period of 568 days. Histogram shows Hann windowed power 
(a = 2) power of spline interpolation flattened assuming an power-law index 
of -3. 



